Spatiotemporal Correlations of Earthquakes 
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Statistical properties of earthquakes are studied both by the analysis of real earthquake catalog of 
Japan and by numerical computer simulations of the spring-block model in both one and two di- 
mensions. Particular attention is paid to the spatiotemporal correlations of earthquakes, e.g., the 
recurrence-time distribution or the time evolution before and after the mainshock of seismic distri- 
bution functions, including the magnitude distribution and the spatial seismic distribution. Cer- 
tain eminent features of the spatiotemporal correlations, e.g., foreshocks, aftershocks, swarms and 
doughnut-like seismic pattern, are discussed. 



1 Introduction 

Although earthquakes are obviously complex phe- 
nomena, the basic physical picture of earthquakes 
seems to have been well established now: Earth- 
quake is a stick-slip frictional instability of a fault 
driven by steady motions of tectonic plates ^ . 
Although it remains to be extremely difficult at 
the present stage to say something really credi- 
ble for each individual earthquake event, if one 
collects many events and take average over these 
events, a clear tendency often shows up there. 
Thus, it is sometimes possible to say something 
credible for the average or statistical properties 
of earthquakes, or more precisely, sets of earth- 
quakes. 

In this article, we wish to review some of 
our recent studies on the statistical properties of 
earthquakes. Our study of earthquakes is moti- 
vated by the following three issues. 

Critical versus characteristic 

It has long been known empirically that cer- 
tain power-laws often appear in the statistical 
properties of earthquakes, e.g., the Gutenberg- 
Richter (GR) law for the magnitude distribution 
of earthquakes, or the Omori law for the time evo- 
lution of the frequency of aftershocks [2]. Power- 
law means that there is no characteristic scale in 
the underlying physical phenomenon. In statis- 
tical physics, one of the most widely recognized 
occasion of the appearance of a power-law is criti- 
cal phenomena associated with a thermodynamic 
second-order phase transition. Indeed, inspired by 
this analogy, Bak and collaborators introduced the 
concept of the "self-organized criticality (SOC)" 
into earthquakes |S1 15 • According to this picture. 
The Earth's crust is always in the critical state 
generated dynamically, and power-laws associated 
with statistical properties of earthquakes are re- 
garded as realizations of the intrinsic critical na- 
ture of earthquakes. Indeed, the SOC idea gives a 



natural explanation of the scale-invariant power- 
law behaviors frequently observed in earthquakes, 
including the GR law and the Omori law. 

In contrast to such an SOC view of earth- 
quakes, an opposite view has also been common 
in earthquake studies, a view which regards earth- 
quakes as "characteristic". In this view, earth- 
quakes are supposed to possess their own charac- 
teristic scales, e.g., a characteristic energy scale or 
a characteristic time scale. 

Thus, whether an earthquake is critical or 
characteristic remains to be one of central issues 
of modern earthquake studies. 

Possible precursory phenomena — spatiotemporal 
correlations of earthquakes 

In conjunction with earthquake prediction, 
possible precursory phenomena associated with 
large earthquakes have special importance. If one 
takes a statistical approach, a natural quantity to 
be examined might be spatiotemporal correlations 
of seismicity, i.e., how earthquakes correlate in 
space and time. If we could identify the property 
in which a clear anomaly is observed preceding the 
large event, it might be useful for earthquake pre- 
diction. In the present article, we wish to investi- 
gate among others various types of spatiotemporal 
correlation functions of earthquakes. 

The constitutive relation and the nature of stick- 
slip dynamics 

Since earthquakes can be regarded as a stick- 
slip frictional instability of a pre-existing fault, 
the statistical properties of earthquakes should be 
governed by the physical law of rock friction |^ 
12]. Unfortunately, our present understanding 
of physics of friction is still poor. We do not 
have precise knowledge of the constitutive relation 
governing the stick-slip dynamics at earthquake 
faults. The difficulty lies partly in the fact that a 
complete microscopic theory of friction is still not 
available, but also in the fact that the length and 



time scales relevant to earthquakes are so large 
that the applicability of laboratory experiments 
of rock friction is not necessarily clear. 

A question of fundamental importance in 
earthquake studies might be how the properties 
of earthquakes depend on the constitutive rela- 
tion and the material parameters characterizing 
earthquake faults, including the elastic properties 
of the crust and the constitutive parameters char- 
acterizing the friction force. 

To answer this question and to get deeper in- 
sight into the physical mechanism governing the 
stick-slip process of earthquakes, a proper model- 
ing of an earthquake might be an important step. 
In fact, earthquake models of various levels of sim- 
plifications have been proposed in geophysics and 
statistical physics, and their statical properties 
have been extensively studied mainly by means 
of numerical computer simulations. One of the 
standard model is the so-called spring-block model 
originally proposed by Burridge and Knopoff (BK 
model) 0, which we will employ in the present 
particle. Yet, our present understanding of the 
question how the earthquake properties depend on 
the constitutive relation and the material param- 
eters characterizing earthquake faults remains far 
from satisfactory. 

In the present article, in order to approach the 
three goals mentioned above, we take two com- 
plementary approaches: In one, we perform nu- 
merical computer simulations based on the spring- 
block model to clarify the spatiotemporal correla- 
tions of seismic events. Both the one-dimensional 
(ID) BK model [HI E] and the two-dimensional 
(2D) BK model are studied. In the other, we 
analyze the earthquake catalog of Japan to ex- 
amine the spatiotemporal correlations of real seis- 
micity [Sj. We then compare the results of nu- 
merical model simulation and the analysis of real 
earthquake catalog, hoping that such a compar- 
ison might give us useful information about the 
nature of earthquakes. 

The following part of the article is organized 
as follows. In section 2, we introduce the model 
employed in our numerical computer simulations, 
and explain some the details of the simulations. 
We also introduce the earthquake catalog used in 
our analysis of real seismicity of Japan. Then, in 
section 3, we report on the results of our analysis 
of the spatiotemporal correlations of real earth- 
quakes based on the earthquake catalog of Japan, 
together with the results of our numerical simula- 
tions of the ID and 2D BK models. The statistical 
properties examined in this section include; (i) the 
magnitude distribution, (ii) the local recurrence- 
time distribution, (iii) the global recurrence-time 



distribution, (iv) the time evolution of the spatial 
seismic distribution before the mainshock, (v) the 
time evolution of the spatial seismic distribution 
after the mainshock, and (vi) the time-resolved 
magnitude distribution before and after the main- 
shock. Finally, section 4 is devoted to summary 
and discussion of our results. 

2 The model, the simulation and 
the catalog 

In this section, we introduce the spring-block 
model which we will use in our numerical com- 
puter simulation, and explain some of the details 
of the simulation. We also introduce the seismic 
catalog of Japan which we will use in our analysis 
of the spatiotemporal correlations of real earth- 
quakes. The results of our numerical computer 
simulations and the analysis of seismic catalog of 
Japan will be presented in the following section 3. 

2.1 The spring-block model of earth- 
quakes 

The spring-block model of earthquakes was orig- 
inally proposed by Burridge and Knopoff [3]. In 
this model, an earthquake fault is simulated by 
an assembly of blocks, each of which is connected 
via the elastic springs to the neighboring blocks 
and to the moving plate. All blocks are subject 
to the friction force, the source of the nonlin- 
earity in the model, which eventually realizes an 
earthquake-like frictional instability. The model 
contains several parameters representing, e.g., the 
elastic properties of the crust and the frictional 
properties of faults. 

In the ID BK model, the equation of motion 
for the i-th block can be written as 

mtii = kp{u't' - Ui) + K{Ui+i - 2Ui + Ui-i) - ^i, 

(1) 

where t' is the time, Ui is the displacement of the 
i-ih. block, v' is the loading rate representing the 
speed of the plate, and <I>j is the friction force 
working at the block i. 

In order to make the equation dimensionless, 
we measure the time t' in u nits of the character- 
istic frequency w = kp/m and the displacement 

Ui in units of the length L = ^o/kp, being the 
static friction. Then, the equation of motion can 
be written in the dimensionless form as 

Ui = ut-Ui + l^{ui+i - 2ui + nj_i) - (j)i, (2) 

where t = t'to is the dimensionless time, Ui = Ui/ L 
is the dimensionless displacement of the z-th block. 



/ = \Jkc/kp is the dimensionless stiffness parame- 
ter, u = v' l(Lijj) is the dimensionless loading rate, 
and (/>i = ^ij^Q is the dimensionless friction force 
working at the block i. 

The form of the friction force (/> is specified by 
the constitutive relation. As mentioned, this part 
is a vitally important, yet largely ambiguous part 
in the proper description of earthquakes. In order 
for the model to exhibit a dynamical instability 
corresponding to an earthquake, it is essential that 
the friction force </> possesses a frictional weakening 
property, i.e., the friction should become weaker 
as the block slides. 

As the simplest form of the friction force, 
we assume here the form used by Carlson et al, 
which represents the velocity-weakening friction 
force jlUl lllj . Namely, the friction force </>(«) is 
assumed to be a single- valued function of the ve- 
locity Ui, i.e., (pi, gets smaller as Uj increases, 

.X / (-00,1], for iii < 0, 

nu) = < 1^ f ^ Q (3) 

where its maximum value corresponding to the 
static friction has been normalized to unity. As 
noted above, this normalization condition (j){u = 
0) = 1 has been utilized to set the length unit L. 
The back-slip is inhibited by imposing an infinitely 
large friction for Uj < 0, i.e., 4){u < 0) = — cx3. 

In this velocity-weakening constitutive rela- 
tion, the friction force is characterized by the two 
parameters, a and a. The former, a, represents 
an instantaneous drop of the friction force at the 
onset of the slip, while the latter, a, represents 
the rate of the friction force getting weaker on in- 
creasing the sliding velocity. In our simulation, we 
regard a to be small, and fix cr = 0.01. 

It should be emphasized again that, although 
the above Carlson-Langer velocity-weakening fric- 
tion is rather simple and has been widely used in 
numerical simulations, the real constitutive rela- 
tion might not be so simple with features possi- 
bly different from the simplest velocity-weakening 
one. Indeed, there have been several other pro- 
posals for the law of rock friction, e.g., the slip- 
weakening friction force ^ 2, 12 , _13, or the rate- 
and state-dependent friction force ^ El 
EIEIE]- In this article, we leave the study of 
these different constitutive relations to other refer- 
ences or to future studies, and assume the simplest 
velocity- weakening friction force given above. 

We also assume the loading rate v to be in- 
finitesimally small, and put v = {) during an earth- 
quake event, a very good approximation for real 
faults jlOl lllj . Taking this limit ensures that the 
interval time during successive earthquake events 



can be measured in units of v irrespective of par- 
ticular values of v. Taking the — > limit also 
ensures that, during an ongoing event, no other 
event takes place at a distant place, independently 
of this ongoing event. 

The extension of the ID BK model to 2D is 
rather straightforward. In 2D, the blocks are con- 
sidered to be arranged in the form of a square ar- 
ray connected with the springs of the spring con- 
stant kc We consider the isotropic and uni- 
form case where kc is uniform everywhere in the 
array independent of the spatial directions. All 
blocks are connected to the moving plate via the 
springs of the spring constant kp, and are also sub- 
ject to the velocity- weakening friction force de- 
fined above. The plate is driven along the x- 
direction with a constant rate v. It is assumed 
that all blocks move along the x-direction only, 
i.e., the displacement of the block at the posi- 
tion {ix,iy) is assumed to be given by u{ix,iy) = 

{Ux{ix,iy),0)- 

Since the early study by Burridge and Knopoff 
1^, the properties of this BK model has been stud- 
ied, simulations, 

Carlson, Langer and collaborators performed 
a pioneering study of the statistical properties of 
the ID BK model quite extensively [Till ITTl [THl 
'20' '21', with particular attention to the magnitude 
distribution of earthquake events. It was observed 
that, while smaller events persistently obeyed the 
GR law, i.e., staying critical or near-critical, larger 
events exhibited a significant deviation from the 
GR law, being off-critical or "characteristic" [TUl 
19, 20 . 211 0^. Shaw, Carlson and Langer studied 
the same model by examining the spatiotemporal 
patterns of seismic events preceding large events, 
observing that the seismic activity accelerates as 
the large event approaches j23j . 

The BK model was also extended in several 
ways, e.g., taking account of the effect of the vis- 
cosity |241 1251 126j , modifying the form of the fric- 
tion force |241 126j , or taking account of the long- 
range interactions PZj . The 2D version of the BK 
model was also analyzed by Carlson and by 
Myers et al 

The study of statistical properties of earth- 
quakes was promoted in early nineties, inspired 
by the work by P. Bak and collaborators who 
emphasized the concept of "self-organized criti- 
cality (SOC)" HI. The SOC idea was devel- 
oped mainly on the basis of the cellular-automaton 
versions of the earthquake model jSl IS I2H1 HSl 
EOl EH El IHI- The statistical properties of 
these cellular-automaton models were also stud- 
ied quite extensively, often interpreted within the 
SOC framework. These models apparently repro- 
duce several fundamental features of earthquakes 



such as the GR law, the Omori law of after- 
shocks, the existence of foreshocks, etc. We note 
that, although many of these cellular-automaton 
models were originally introduced to mimic the 
spring-block model, their statistical properties 
are not always identical with the original spring- 
block model. Furthermore, as compared with the 
spring-block model, the cellular-automaton mod- 
els are much more simplified so that the model 
does not have enough room to represent vari- 
ous material properties of the earthquake fault 
in a physically appealing way. Thus, the spring- 
block model has an advantage over the cellular- 
automaton models that the dependence on the 
material parameters, including the constitutive 
and elastic properties, are more explicit. 

2.2 The numerical simulation 

As already mentioned, numerical model simula- 
tion is a quite useful tool for our purposes: First, 
generating huge number of events required to at- 
tain the high precision for discussing the statistical 
properties of earthquakes is easy to achieve in nu- 
merical simulations whereas it is often difficult to 
achieve in real earthquakes, especially for larger 
ones. Second, various material parameters char- 
acterizing earthquake faults are extremely difficult 
to control in real faults, whereas they are easy to 
control in model simulations. 

We solve the equation of motion (2) by us- 
ing the Runge-Kutta method of the fourth order. 
The width of the time discretization At is taken 
to be Atu = 10~^. We have checked that the 
statistical properties given below are unchanged 
even if we take the smaller At. Total number of 
10^ events are generated in each run, which are 
used to perform various averagings. In calculating 
the observables, initial 10^ events are discarded as 
transients. 

In order to eliminate the possible finite-size ef- 
fects, the total number of blocks N are taken to be 
large. In ID, we set N = 800, imposing the peri- 
odic boundary condition. The size dependence of 
the results was examined in Ref . |7j with varying 
in the range 800 < N < 6400. In 2D, we set our 
lattice size 160 x 80 with the periodic boundary 
condition on the longer side, and the free bound- 
ary condition on the shorter side. 

We study the properties of the model, with 
varying the frictional parameter a and the elas- 
tic parameter I. In this article, attention is paid 
to the dependence on the parameter a, since the 
parameter a, which represents the extent of the 
frictional weakening, turns out to affect the result 
most significantly [7]. 



2.3 The seismic catalog of Japan 

In order to compare the simulation data on 
the ID and 2D BK models with the corre- 
sponding data for real earthquakes, we analyze 
in the following section the seismic catalog 
of Japan provided by Japan University Net- 
work Earthquake Catalog (JUNEC) available at 



http:/ /kea.eri.u-tokyo. ac.jp/CATALOG/junec/monthly.html 
The catalog covers earthquakes which occurred 
in Japan area during July 1985 and December 
1998. Total of 199,446 events are recorded in 
the catalog. As an example, we show in Fig^ 
a seismicity map of Japan generated from the 
JUNEC catalog, where all large earthquakes 
of their magnitudes greater than five, which 
occurred in Japan area during 1985-1998, are 
mapped out 
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Figure 1: A seismicity map of Japan generated from 
the JUNEC catalog. Large earthquakes of their magni- 
tudes greater than five, which occurred in Japan area 
during 1985-1998, are mapped out by using the pro- 
gram^Jleveloped by H. Tsuruoka [55]. 



3 Statistical properties of earth- 
quakes 

In this section, we show the results of our numeri- 
cal simulations of the ID and 2D BK models, and 
compare them with the results of our analysis of 



the seismic catalog of Japan (JUNEC catalog). 
We study several observables, i.e., (i) the mag- 
nitude distribution, (ii) the local recurrence-time 
distribution, (iii) the global recurrence-time dis- 
tribution, (iv) the time evolution of the spatial 
seismic distribution before the mainshock, (v) the 
time evolution of the spatial seismic distribution 
after the mainshock, and (vi) the time-resolved 
magnitude distribution before and after the main- 
shock. We show these results consecutively below. 
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Figure 2: The magnitude distribution of earthquakes 
in Japan generated from the JUNEC catalog. The data 
for m > 3 lie on a straight line with the GR exponent 
&~0.9. 



3.1 The magnitude distribution 

In Fig|21 we show the magnitude distribution 
R{m) of earthquakes of Japan generated from the 
JUNEC catalog, where R{m)dm represents the 
rate of events with their magnitudes in the range 
[m, m +dm]. The data lie on a straight line fairly 
well for m > 3, obeying the GR law. The slope 
of the straight line gives the power-law exponent 
about b ~ 0.9. 

In FiglHfa), we show the magnitude distribu- 
tion of earthquakes calculated from our nu- 
merical simulation of the ID BK model. The pa- 
rameter a is varied in the range 0.25 < a < 10, 
while the elastic parameter I is fixed to / = 3. 

In the BK model, the magnitude of an event, 
^, is defined as the logarithm of the moment Mq, 
i.e., 

^l = lnMo, Mo = ^Aui, (4) 

i 

where Auj is the total displacement of the i-th 
block during a given event and the sum is taken 
over all blocks involved in the event [10 . It should 
be noticed that the absolute numerics of the mag- 
nitude value of the BK model fj, has no direct 
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Figure 3: The magnitude distribution of seismic 
events calculated for the ID BK model, for the range 
of larger a > 1 (a), and for the range of smaller a < 1 
(b). 

quantitative correspondence to the magnitude m 
of real earthquakes. 

As can be seen from FigOl the form of the 
calculated magnitude distribution depends on the 
a-value considerably. The data for a = 1 lie on a 
straight line fairly well, apparently satisfying the 
GR law. The value of the exponent B describing 
the GR-like power-law behavior, oc 10~^, is esti- 
mated to be ~ 0.50 corresponding to 6 ~ 0.75, 
which is slightly smaller than the 6- value observed 
for the JUNEC catalog b ~ 0.9. Remember the 
relation b = ^B. 

On the other hand, the data for larger a, i.e., 
the ones for a > 2 deviate from the GR law at 
larger magnitudes, exhibiting a pronounced peak 
structure, while the power-law feature still re- 
mains for smaller magnitudes. These features of 
the magnitude distribution are consistent with the 
earlier observation of Carlson and Langer jlOl I20j . 
It means that, while smaller events exhibit self- 
similar critical properties, larger events tend to ex- 
hibit off-critical or characteristic properties, much 
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Figure 4: The magnitude distribution of seismic 
events calculated for the 2D BK model. 

more so as the velocity-weakening tendency of the 
friction is increased. The observed peak structure 
gives us a criterion to distinguish large and small 
events. Below, we regard events with their mag- 
nitudes fi greater than /Xc = 3 as large events of 
the ID BK model, /ic = 3 being close to the peak 
position of the magnitude distribution of Fig|2[a). 
In an earthquake with ^ = 3, the mean number 
of moving blocks are about 76 (a = 1) and 60 
(q = 2,3). 

As can be seen from Fig|2Ib), the data at 
smaller a < 1 exhibit considerably different be- 
haviors from those for a > 1. Large events are 
suppressed here. For a = 0.25, in particular, all 
events consist almost exclusively of small events 
only. This result might be consistent with the ear- 
lier observation which suggested that the smaller 
value of Q < 1 tended to cause a creeping-like be- 
havior without a large event |2n|. In particular, 
Vasconcelos showed that a single block system ex- 
hibited a "first-order transition" at a = 0.5 from 
a stick-slip to a creep , whereas this discontin- 
uous transition becomes apparently continuous in 
many-block system |371 138| . Since we are mostly 
interested in large seismic events here, we con- 
centrate in the following on the parameter range 
a > 1. In contrast to the parameter a, the mag- 
nitude distribution turns out to be less sensitive 
to the stiffness parameter I. Further details of the 
/-dependence is given in Ref.0. 

In Fig2J we show the magnitude distribution 
R{^) of the 2D BK model with varying the a 
value. At larger magnitudes, a deviation from 
the GR law is observed for any value of a, i.e., 
the calculated magnitude distribution exhibits a 
peak structure at larger magnitude irrespective 
of the a- value, suggesting that larger earthquakes 
tend to be characteristic. From the observed peak 



structure, we regard events with their magnitudes 
/i greater than //c = 5 as large events of the 2D 
BK model. 

While the BK model tends to reproduce the 
GR law for earthquakes of smaller magnitudes, its 
relevance to the GR law observed in real seimicity 
is not entirely clear. It has been suggested that 
the GR 6-value might be related to the fractal di- 
mension of the fault interface |391 1401 141j . More 
recently, Chakrabarti and collaborators proposed 
an earthquake model in which the magnitude dis- 
tribution of earthquakes is related to the contact 
area distribution between the two fractal surfaces 
of the plates [l^l- In these approaches, the in- 
trinsic nonuniformity at earthquake faults plays 
an essential role in realizing the GR law and the 
critical nature of earthquakes. Whether this is re- 
ally true, as well as its possible relation to the BK 
model where the nonuniformity is apparently ab- 
sent, or at least not explicit, needs to be examined 
further. 

3.2 The local recurrence-time distribu- 
tion 

A question of general interest may be how large 
earthquakes repeat in time, do they occur near pe- 
riodically or irregularly? One may ask this ques- 
tion either locally, i.e., for a given finite area on 
the fault, or globally, i.e., for an entire fault sys- 
tem. The picture of characteristic earthquake pre- 
sumes the existence of a characteristic recurrence 
time. In this case, the distribution of the recur- 
rence time of large earthquakes, T, is expected to 
exhibit a peak structure at such a characteristic 
time scale. If the SOC concept applies to large 
earthquakes, by contrast, such a peak structure 
would not show up. 

In the upper panel of Fig|Sl we show the dis- 
tribution of the local recurrence time T of large 
earthquakes of Japan with m > rric = 4, calcu- 
lated from the JUNEC catalog. In defining the 
recurrence time locally, the subsequent large event 
is counted when a large event with m > rric = 4 
occurs with its epicenter lying within a circle of 
radius 30 km centered at the epicenter of the pre- 
vious large event. The mean recurrence time T is 
then estimated to be 148 days. In the middle and 
lower panels of Figs El the same data are re-plotted 
on a double-logarithmic scale [middle panel], and 
on a semi- logarithmic scale [lower panel] , with the 
magnitude threshold rric = 3,4 and 5. One sees 
from these figures that the local recurrence-time 
distribution exhibits power-law-like critical fea- 
tures at shorter times, which seems to cross over 
to a faster exponential-like decay at longer times. 
The time range in which the data obey a power- 
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Figure 5: The local recurrence-time distribution of 
large earthquakes in Japan generated from the JUNEC 
catalog. The distribution is given on a bare scale for 
large earthquakes with rric = 4 [upper figure], on a 
logarithmic scale with rUc = 3,4 and 5 [middle figure], 
and on a semi- logarithmic scale with nic = 3,4 and 5 
[lower figure]. 



law tends to get longer for larger earthquakes. The 
power-law exponent estimated for earthquakes of 
rric = 5 in the time range T < 1000 days is about 
|. This value is a bit smaller than the standard 
Omori exponent ~ 1. 

In the upper panel of Figl^l we show the dis- 
tribution of the local recurrence time T of large 
earthquakes with n > fJ-c = 'i, calculated for the 
ID BK model. In the insets, the same data in- 
cluding the tail part are re-plotted on a semi- 
logarithmic scale. In defining the recurrence time 
locally in the BK model, the subsequent large 
event is counted when a large event occurs with 
its epicenter in the region within 30 blocks from 
the epicenter of the previous large event. The 
mean recurrence time T is then estimated to be 
fu = 1.47, 1.12, and 1.13 for a = 1, 2 and 3, 
respectively. 

As can be seen from FigEl the tail of the dis- 
tribution is exponential at longer T > T for all 
values of a, while the form of the distribution at 
shorter T < T is non-exponential and differs be- 
tween for a = 1 and for a = 2 and 3. For a = 2 
and 3, the distribution has an eminent peak at 
around Tv ~ 0.5, not far from the mean recur- 
rence time. This means the existence of a char- 
acteristic recurrence time, suggesting the near- 
periodic recurrence of large events. This charac- 
teristic behavior is in sharp contrast to the critical 
behavior without any peak structure observed in 
the JUNEC data. It should also be mentioned, 
however, that there were reports in the literature 
of a near-periodic recurrence of large events at sev- 
eral real faults PUS]. 

For a = 1, by contrast, the peak located close 
to the mean T is hardly discernible. Instead, the 
distribution has a pronounced peak at a shorter 
time Ti/ ~ 0.10, just after the previous large event. 
In other words, large events for a = 1 tend to oc- 
cur as "twins". This has also been confirmed by 
our analysis of the time record of large events. In 
fact, as shown in Ref.[7], a large event for the case 
of a = 1 often occurs as a "unilateral earthquake" 
where the rupture propagates only in one direc- 
tion, hardly propagating in the other direction. 
When a large earthquake occurs in the form of 
such a unilateral earthquake, further loading due 
to the plate motion tends to trigger the subse- 
quent large event in the opposite direction, caus- 
ing a twin-like event. This naturally explains the 
small-T peak observed in FigslHlfor a = 1. 

In the lower panel of FigEl the local 
recurrence-time distribution of large events is 
shown for the cases of q = 1, with varying the 
magnitude threshold as /Uc = 2, 3 and 4. As can 
be seen form this figure, the form of the distribu- 
tion for a = 1 largely changes with the threshold 




T/T 

Figure 6: The local recurrence-time distribution of 
large events of the ID BK model. In the upper figure, 
the distribution is given for a = 1,2 and 3 with fixing 
lie = 3, whereas, in the lower figure, the distribution 
is given with varying the magnitude threshold as fic — 

2, 3 and 4 with fixing a = 1. 

value /ic- Interestingly, in the case of /ic = 4, 
the distribution has two distinct peaks, one corre- 
sponding to the twin-like event and the other to 
the near-periodic event. Thus, even in the case 
of a = 1 where the critical features are appar- 
ently dominant for smaller thresholds /ic = 2 and 

3, features of a characteristic earthquake becomes 
increasingly eminent when one looks at very large 
events. 

In Fig|7J the local recurrence-time distribution 
is shown for the 2D BK model. In defining the re- 
currence time locally in the 2D BK model, a sub- 
sequent large event is counted when a large event 
occurs with its epicenter lying within a circle of 
its radius 5 blocks centered at the epicenter of the 
previous large event. The mean recurrence time T 
is then estimated to be Tv = 1.47, 1.12, and 1.13 
for a = 1, 2 and 3, respectively. The behavior 
of the 2D model is similar to the ID model, with 
enhanced characteristic features. The peak struc- 
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Figure 7: The local recurrence-time distribution of 
large events of /i > /ic = 5 of the 2D BK model. 

ture of the distribution is very prominent even for 
a = 1. 

3.3 The global recurrence-time distri- 
bution 

In this subsection, we examine the global 
recurrence-time distribution. By the word 
"global", we consider the situation where the re- 
gion in which one identifies the next event is suf- 
ficiently wide, much wider than the typical size of 
the rupture zone of large events. 

From the JUNEC catalog, we construct in the 
upper panel of FigEJ such global recurrence-time 
distribution for entire Japan for large earthquakes 
with nic = 4. The mean recurrence time T is 
then estimated to be 0.73 days. In the middle and 
lower panels of Figs El the same data are re-plotted 
on a double-logarithmic scale [middle], and on a 
semi- logarithmic scale [lower] , with the magnitude 
threshold nic = 3,4 and 5. One sees from these fig- 
ures that the global recurrence-time distribution 
of the JUNEC data is very much similar to the 
corresponding local recurrence-time distribution 
shown in FigsEl It exhibits power-law-like critical 
features at shorter times, which crosses over to a 
faster exponential-like decay at longer times. The 
time range in which the data obey a power-law 
gets longer for larger earthquakes. The exponent 
describing the power-law regime is roughly about 
|, which is not far from the corresponding expo- 
nent of the local recurrence-time distribution. 

FigEl exhibits the global recurrence-time dis- 
tribution of large events with //c = 3 calculated 
for the ID BK model. As can clearly be seen from 
the figure, in the BK model, the form of the dis- 
tribution takes a different form from the local one: 
The peak structure seen in the local distribution 
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Figure 8: The global recurrence-time distribution of 
large earthquakes in entire Japan generated from the 
JUNEC catalog. The distribution is given on a bare 
scale for large earthquakes with rric = 4 [upper figure] , 
on a double-logarithmic scale with varying rric = 3,4 
and 5 [middle figure], and on a semi- logarithmic scale 
with varying rric = 3,4 and 5 [lower figure]. 



no longer exists here. Furthermore, the form of 
the distribution tail at larger T is not a simple 
exponential, faster than exponential: See a cur- 
vature of the data in the inset of FigEl These 
features of the global recurrence-time distribution 
turn out to be rather robust against the change of 
the parameter values such as a and I, as long as 
the system size is taken to be sufficiently large. 

The observation that the local and the global 
recurrence-time distributions exhibit mutually dif- 
ferent behaviors means that the form of the dis- 
tribution depends on the length scale of mea- 
surements. Such scale-dependent features of the 
recurrence-time distribution of the BK model 
are in apparent contrast with the scale-invariant 
power-law features of the recurrence-time distri- 
butions observed in the JUNEC data shown in 
FigsEl and |S1 and the ones reported for some of 
real faults [441 145j . We note that essentially the 
same behavior as the one of the ID BK model is 
also observed in the 2D BK model (the data not 
shown here). 
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Figure 9: The global recurrence-time distribution of 
large events of the ID BK model. The distribution 
is given for a — 1,2 and 3 with fixing the magnitude 
threshold /ic — 3. 



3.4 Spatiotemporal seismic correla- 
tions before the mainshock 

In this subsection, we study the spatiotemporal 
correlations of earthquake events before the main- 
shock. 

We begin with the JUNEC catalog. First, we 
define somewhat arbitrarily the mainshock as an 
event whose magnitude m is greater than rUc = 5, 
and pay attention to the frequency of earthquake 
events preceding the mainshock, particularly its 
time and distance dependence. Both the time t 



and the distance r are measured relative to the 
time and the position of the subsequent main- 
shock. The distance between the two events is 
measured here as the distance between their epi- 
centers. The event frequency is normahzed by the 
factor r associated with the area element of the 
polar coordinate. 
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Figure 10: The time evolution of the event frequency 
of arbitrary magnitude preceding the mainshock, plot- 
ted as a function of the distance from the epicenter 
of the upcoming mainshock. The magnitude threshold 
for the mainshock is = 5 (a) and — 3 (b). The 
data are generated from the JUNEC catalog, while the 
contribution from the two special regions, i.e., the Izu 
region and the Ebino region, has been omitted from 
the data: See the text for further details. 

FigllOl exhibits the event frequency plotted as 
a function of r, the distance from the upcom- 
ing mainshock, for four time periods preceding 
the mainshock each containing five days. The 
data show the time evolution of the spatial pat- 
tern of event number irrespective of its magni- 
tude, which occur within 30 km from the epicenter 
of the upcoming mainshock and during the last 
20 days toward the mainshock. The data have 
been averaged over the mainshocks contained in 



the data set taken over entire Japan, but omitting 
the contributions of the two special narrow regions, 
i.e., Izu and Ebino. In these two regions, active 
earthquake "swarms" occurred during the period, 
which lead to significantly different behavior in 
the spatiotemporal correlations (we return to this 
point later in this subsection). In the data set, 
total of 990 mainshocks are included, each main- 
shock accompanying about 8.7 preceding events 
on average in the time/space range shown in the 
figure. As can clearly be seen from Fig llOf a). there 
is a tendency that the seismic activity accelerates 
as the mainshock approaches, and this tendency is 
more pronounced in a closer vicinity of the epicen- 
ter of the upcoming mainshock. The data demon- 
strate that the mainshock accompanies foreshocks, 
at least in the statistical sense. 

It should be emphasized here that the evidence 
of foreshock activity becomes clear only after av- 
eraging over a large number of mainshocks. Note 
that, in the time/space region of interest, each sin- 
gle mainshock accompanies on average less than 
ten foreshocks only, a too small number to say 
anything definite. 

One may wonder if certain qualitative fea- 
tures of the spatiotemporal correlations shown 
in Fig llOf a) might change depending on the re- 
gion, the depth and the magnitude-threshold of 
the mainshock. Examination of the JUNEC data, 
however, reveals that the qualitative features of 
Fig llOf a) are rather robust against these param- 
eters. The average number of foreshocks depends 
somewhat on each region in Japan: For example, 
the mean numbers of foreshocks in the space /time 
region of Fig llOf a') are 9.4 for the northern part 
of Japan (40° N- latitude), 14.7 for the middle 
part (36°N-40°N) and 3.6 for the southern part (- 
36°N). Nevertheless, the qualitative features of the 
spatiotemporal correlations turn out to be more or 
less common. 

In Fig llOf b). we show the similar spatiotempo- 
ral seismic correlations as in Fig JlOf a) , but now re- 
ducing the magnitude threshold of the mainshock 
from rric = 5 to rUc = 3. By this definition, the to- 
tal number of mainshocks is increased to 53,835 so 
that the better statistics is expected. Indeed, the 
data shown in Fig llOf b) are far less erratic than 
those in Fig JlOf a). refiecting the improvement of 
the statistics. Nevertheless, the qualitative fea- 
tures remain almost the same as in Fig JlOr a). 

In order to clarify the relation with the inverse 
Omori law which describes the time evolution of 
the frequency of foreshocks, we show in Fig llll the 
time dependence of the frequency of foreshocks be- 
fore the mainshock on a double-logarithmic scale: 
In the upper panel, the data with the distance 
range r^ax = 30 km are shown with varying the 



magnitude threshold as nic = 3,4 and 5, while in 
the lower panel, the data with rUc = 5 are shown 
with varying the distance range as rmax = 5, 30 
and 300 km. Except for the case of very large 
distance range rmax = 300 km, the inverse Omori 
exponent comes around 0.5. 

Now, we wish to turn to the spatiotemporal 
correlations in the two special regions, the contri- 
bution of which have intentionally been omitted 
in our analysis of Figs llOl and 1111 In fact, these 
special regions are associated with "swarms". If 
there occurs an earthquake swarm which contains 
in itself a few large events which can be regarded 
as mainshocks, they make a significant contri- 
bution to the above spatiotemporal correlations 
because a huge number of small events are con- 
tained in swarms. The two swarms we discarded 
in our analysis of Figs llOl and ^2 ^-I'e, (i) the Izu 
earthquake swarm, and (ii) the Ebino earthquake 
swarm (Ebino is located close to the Kirishima vol- 
canoes in southern Kyusyu). Indeed, if the contri- 
bution of these two swarm regions were not sepa- 
rated in our analysis of Figs llOl and llll the result- 
ing distribution function would look considerably 
different. 

In Figs ll2l we show the the spatiotemporal 
seismic correlations associated with mainshocks 
with rric = 5 in the Izu swarm region [34. 8° N- 
35.1°N latitude; 139.0°E-139.5°E longitude] (a), 
and in the Ebino swarm region [31.8°N-32.0°N; 
130.2°E-130.6°E] (b), respectively. The total 
numbers of mainshocks are 4 (Izu) and 6 (Ebino) 
here. In either case, the average number of fore- 
shocks per mainshock is very large. In the case of 
Izu, it is about 256 within the range of 30km and 
20 days from the mainshock, which is an order of 
magnitude larger than the number associated with 
swarm-unrelated mainshocks. 

As can immediately be seen from Figs ll2l the 
spatiotemporal pattern of these regions are quite 
different from the one in other regions shown 
in Figs llOl which suggests that the property of 
swarm-related earthquakes might differ qualita- 
tively from that of swarm- unrelated earthquakes. 
For example, foreshocks in the Izu swarm region 
suddenly accelerate on the onset time of about a 
week, with the center of activity about 5 km away 
from the mainshock position, the associated spa- 
tial correlation function exhibiting a pronounced 
peak at about r ~ 5 km. This seismic pattern 
is not dissimilar to the doughnut-like seismic pat- 
tern discussed by Mogi as occurring preceding the 
mainshock (Mogi doughnut) |^. However, since 
the doughnut-like seismic pattern as observed in 
Figs ll2l is not observed in our analysis of more 
generic case of the wider region, at least in a statis- 
tically significant manner, it is likely to be related 
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Figure 11: The frequency of seismic events of ar- 
bitrary magnitude preceding the mainshock, plotted 
versus the time until the mainshock on a double- 
logarithmic scale. The data are generated from the 
JUNEC catalog, where the contribution from the two 
special regions, i.e., the Izu region and the Ebino re- 
gion, has been omitted: See the text for further details. 
In the upper figure, the magnitude threshold of the 
mainshock is varied as nic = 3, 4 and 5 with fixing the 
distance range of observation Vmax = 30 km, whereas, 
in the lower figure, the distance range of observation 
is varied as rmax = 5, 30 and 300 km with fixing the 
magnitude threshold rric — 5. 
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Figure 12: The time evolution of the event frequency 
of arbitrary magnitude preceding the mainshock, plot- 
ted versus the distance from the epicenter of the up- 
coming mainshock. The data are generated from the 
JUNEC catalog, for the Izu earthquake swarm re- 
gion [34.8°N-35.1°N latitude; 139.0°E-139.5°E longi- 
tude] (a), and for the Ebino earthquake swarm region 
[31.8°N-32.0°N; 130.2°E-130.6°E] (b), both during the 
period 1985-1998. The magnitude threshold of the 
mainshock is m > mr = 5 



to the specialty of Izu (and Ebino) region. Swarm 
activity is considered to be related to the activity 
of underground magma or water. Hence, some of 
the statistical properties of swarm-related earth- 
quakes might be better discriminated from those 
of more general swarm-unrelated earthquakes. 

We examine next the corresponding spa- 
tiotemporal correlations of the BK model. In the 
BK model, the distance from the epicenter is mea- 
sured in units of blocks. Fig J13l reDresents the time 
evolution of the spatial distribution of seismicity 
before the mainshock with /i > /Xc = 3 for a = 1. 
Very much similar behaviors are observed for other 
values of a. Similarly to the JUNEC data, pre- 
ceding the mainshock, there is a tendency of the 
frequency of smaller events to be enhanced at and 



around the epicenter of the upcoming mainshock. 
This was also observed previously by Shaw et al 
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Figure 13: The time evolution of the spatial seismic 
distribution function before the mainshock of /x > /Zc = 
3, calculated for the ID BK model. The inset repre- 
sents a similar plot with longer time intervals. 

Interestingly, however, as the mainshock be- 
comes imminent, the frequency of smaller events 
is suppressed in a close vicinity of the epicenter of 
the upcoming mainshock, though it continues to 
be enhanced in the surroundings (Mogi doughnut) 
IS| • We note that the quiescence observed here 
occurs only in a close vicinity of the epicenter of 
the mainshock, within one or two blocks from the 
epicenter, and only at a time close to the main- 
shock. The time scale for the appearance of the 
doughnuts-like quiescence depends on the o-value: 
The time scale of the onset of the doughnut-like 
quiescence tends to be longer for larger a. It is not 
clear at the present stage whether the doughnut- 
like quiescence as observed here in the BK model 
has any relevance to the one observed in Figs ll2l 
for certain earthquake swarms. 

As was discussed in detail in Ref . j7j , in the BK 
model, the size of the "hole" of the doughnut-like 
quiescence as well as its onset time scale have no 
correlation with the magnitude of the upcoming 
event. In other words, the doughnut-like quies- 
cence is not peculiar to large events in the present 
model. This means that, by monitoring the on- 
set of the "hole" in the seismic pattern of the BK 
model, one can certainly deduce the time and the 
position of the upcoming event, but unfortunately, 
cannot tell about its magnitude. Yet, one might 
get some information about the magnitude of the 
upcoming event, not from the size and the onset 
time of the "hole" , but from the size of the "ring" 
surrounding the "hole". Thus, we show in Figs ll4l 
the spatial correlation functions before the main- 



shock in the time range < tiy < 0.001 for the 
case of a = 2, with varying the magnitude range 
of the upcoming event. In the figure, the direction 
in which the rupture propagates farther in the up- 
coming event is always taken to be the positive 
direction r > 0, whereas the direction in which 
the rupture propagates less is taken to be the neg- 
ative direction r < 0. As can be seen from the 
figure, although the size of the "hole" around the 
origin r = has no correlation with the magni- 
tude of the upcoming event as mentioned above, 
the size of the region of the active seismicity sur- 
rounding this "hole" is well correlated with the 
size and the direction of the rupture of the up- 
coming event. This coincidence might enable one 
to deduce the position and the size of the upcom- 
ing event by monitoring the pattern of foreshocks, 
although it is still difficult to give a pinpoint pre- 
diction of the time of the upcoming mainshock. 
We note that such a correlation between the size of 
the seismically active region and the magnitude of 
the upcoming event was observed in the BK model 
in Ref.|47j, and was examined in real earthquake 
data as well (15] . 



measure factor r associated with the area element 
of the polar coordinate. The time scale of the on- 
set of the doughnut-like quiescence seems to be 
longer in 2D than in ID. 
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Figure 15: The time evolution of the spatial seismic 
distribution function before the mainshock of ii > fic = 
5, calculated for the 2D BK model. The inset repre- 
sents a similar plot with longer time intervals. 
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Figure 14: The frequency of seismic events preceding 
the events of various magnitude range plotted versus r, 
the distance from the epicenter of the upcoming main- 
shock. The curves correspond to the large events of 
^ > fic = 2,3 and 4, and to the smaller events of 
fi < —2. The direction in which the rupture prop- 
agates farther in the upcoming mainshock is always 
taken to be the positive direction r > 0, whereas the 
direction in which the rupture propagates less is taken 
to be the negative direction r < 0. The parameters 
are taken to be a = 2 and I = 3, the time range being 
< ti/ < 0.001 before the mainshock. 

As shown in Fig ll51 a very much similar be- 
havior including the doughnut-like quiescence is 
observed also in the 2D BK model. Note that the 
event frequency has been normalized here by the 



3.5 Spatiotemporal seismic correla- 
tions after the mainshock 

In this subsection, we examine the spatiotemporal 
correlations of earthquake events after the main- 
shock. The time evolution of the spatial seismic 
pattern calculated from the JUNEC catalog are 
shown in Fig |16l with the magnitude threshold for 
the mainshock nic = 5 (a), and nic = 3 (b). As 
can be seen from these figures, aftershock activity 
is clearly observed. The rate of the aftershock ac- 
tivity is highest just at the epicenter of the main- 
shock, not in the surroundings. It is sometimes 
mentioned in the literature that the aftershock ac- 
tivity is highest near the edge of the rupture zone 
of the mainshock (see, e.g., Ref. |2j). However, this 
is not the case here. 

In order to further clarify the relation with the 
Omori law, we show in Fig |17l the time dependence 
of the frequency of aftershocks after the mainshock 
on a double-logarithmic scale: In the upper panel, 
the data with the distance range rmax = 30 km are 
shown with varying the magnitude threshold as 
rric = 3,4 and 5, while, in the lower panel, the data 
with rric = 5 are shown with varying the distance 
range as r^nax = 5, 30 and 300 km. In all cases 
analyzed, a straight-line behavior corresponding 
to the power-law decay is observed (Omori law). 
In contrast to the corresponding quantity before 
the mainshock, the slope of the straight line, i.e.. 




Figure 16: The time evolution of the event frequency 
of arbitrary magnitude after the mainshock, plotted 
versus the distance from the epicenter of the main- 
shock. The magnitude threshold for the mainshock is 
TOc = 5 (a), and uric — 3 (b). The data are generated 
from the JUNEC catalog. 

the value of the Omori exponent, seems to depend 
on the distance range Tjykix tind the magnitude of 
the mainshock rric. As the distance range rmax 
gets smaller and the magnitude of the mainshock 
rric gets larger, the Omori exponent tends to get 
larger. If one looks at the range r^ax = 5 km 
for mainshocks greater than rric = 5, the Omori 
exponent is about 0.71. 

The corresponding spatiotemporal correla- 
tions after the mainshock are calculated for the 
ID BK model, and the results are shown in Fig |18l 
for the cases of a = 1. As can be seen from the 
figure, aftershock activity is clearly observed with 
the maximum rate occurring at the epicenter of 
the mainshock. In contrast to the JUNEC data, 
the r-dependence of the spatial seismic distribu- 
tion is non-monotonic, aftershock activity being 
suppressed in the surrounding region where the 
rupture was largest in the mainshock. The seis- 
micity near the epicenter is kept almost constant 



Figure 17: The frequency of seismic events of arbi- 
trary magnitude after the mainshock, plotted versus 
the time after the mainshock on a double-logarithmic 
scale. The data are generated from the JUNEC cat- 
alog. In the upper figure, the magnitude threshold of 
the mainshock is varied as rric = 3, 4 and 5 with fix- 
ing the distance range of observation rmax = 30 km, 
whereas, in the lower figure, the distance range of ob- 
servation is varied as r^nax = 5, 30 and 300 km with 
fixing the magnitude threshold rric — 5. 

in time for some period after the mainshock, say, 
in the time range tv < 0.03, which is in appar- 
ent contrast to the power-law decay as expected 
from the Omori law: See the insets of Fig ll8l At 
longer time scales, the seismicity near the epicen- 
ter seems to decay gradually, although the decay 
observed here is not a power-law decay. A very 
much similar behavior is observed also for the 2D 
BK model: See Fig ]19l Hence, aftershocks obey- 
ing the Omori law is not realized in the BK model, 
as already reported j20i . This is in apparent con- 
trast to the observation for real faults. 

Such an absence of aftershocks in the BK 
model might give a hint to the physical origin of 
aftershocks obeying the Omori law, e.g., they may 
be driven by the slow chemical process at the fault. 
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Figure 18: The time evolution of the spatial seismic 
distribution function after the mainshock of > = 
3, calculated for the ID BK model. The inset repre- 
sents a similar plot with longer time intervals. 
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Figure 20: The time-resolved local magnitude distri- 
bution of earthquakes in Japan before the mainshock 
of m > rric = 4, generated from the JUNEC catalog. 
Preceding the mainshock, the 6- value decreases consid- 
erably from the long-time average value. 

or by the elastoplaciticity associated with the as- 
cenosphere, etc, which are not taken into account 
in the BK model. 
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Figure 19: The time evolution of the spatial seismic 
distribution function after the mainshock of > /ic = 
5, calculated for the 2D BK model. The inset repre- 
sents a similar plot with longer time intervals. 



3.6 The time-dependent magnitude 
distribution 

As an other signature of precursory phenomena, 
we examine a "time-resolved" local magnitude dis- 
tribution for several time periods before the large 
event. Fig |2UI represents such a local magnitude 
distribution before large earthquakes, calculated 
from the JUNEC catalog, i.e., the distribution of 
seismicity in the vicinity rmax = 30 km of the 
epicenter of the mainshock with its magnitude 
greater than rric = 4. It can clearly be seen from 
the figure that the GR law persists even just be- 
fore the mainshock, whereas, preceding the main- 
shock, the GR exponent b gets smaller compared 
with the space- and time-averaged 6-value. For ex- 
ample, 100 days before the mainshock, the 6- value 
becomes about 0.60, considerably smaller than the 
averaged value b ~ 0.88. Such a decrease of the 
6- value was also reported in the literature |^ ISU] . 

For comparison, we show the corresponding lo- 
cal magnitude distributions before the mainshock 
of the BK models both in ID (FigEIl) and in 2D 
fFig |22() . The parameter a is taken to be a = 1. 
In ID, only events with their epicenters within 30 
blocks from the epicenter of the upcoming main- 
shock are counted, while, in 2D, only events with 
their epicenters lying in a circle of its radius 5 
blocks centered at the epicenter of the upcoming 
mainshock are counted. As can be seen from the 
figures, as the mainshock approaches, the form of 
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regime slightly increases first, and then, increases 
just before the mainshock. 
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Figure 21: The time-resolved local magnitude distri- 
bution of the ID BK model before the mainshock of 
/I > /ic = 3. 



\a=l 



1=3 



-4 



3. 

c 



rv=0~0.01 - 
fv=0.01~0.1 - 
fv=0.1~l . 
tv=l~2- 
all / 



-10 



A2 ■ 



-14 




-6 



-4 



-2 



2 



10 



Figure 22: The time-resolved local magnitude distri- 
bution of the 2D BK model before the mainshock of 
H> He ^ 5. 

the magnitude distribution changes significantly. 
In ID, the apparent i?-value describing the power- 
law regime tends to increase as the mainshock ap- 
proaches, from the time-averaged value B ~ 0.50 
{b ~ 0.75) to the value B ~ 1.0 {b ~ 1.5) just 
before the mainshock: It is almost doubled. This 
tendency is opposite to what we have just found 
for the JUNEC catalog and several other obser- 
vations for real faults |491 150j . However, a sim- 
ilar increase of the apparent S-value preceding 
the mainshock was reported for some of real faults 
[HT] . For the case of larger a > 1 (the data not 
shown here) , the change of the -B-value preceding 
the mainshock is still appreciable, though in a less 
pronounced manner. More complicated behavior 
is observed in 2D. As the mainshock approaches, 
the apparent -B-value describing the power-law 



Figure 23: The time-resolved local magnitude distri- 
bution of earthquakes in Japan after the mainshock of 
m > rric — 4, generated from the JUNEC catalog. 

Next, we analyze similar time-resolved local 
magnitude distributions, but after the large event. 
Fig |23l represents such a local magnitude distri- 
bution after the large event calculated from the 
JUNEC catalog. In this case, the deviation from 
the averaged distribution is relatively small as 
compared with the one observed before the main- 
shock, although there seems to be a tendency that 
the observed 6-value decreases slightly after the 
mainshock. 
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Figure 24: The time-resolved local magnitude distri- 
bution of the ID BK model after the mainshock of 
/i > /ic = 3. 

For comparison, we show the corresponding lo- 
cal magnitude distributions after the mainshock 
for the BK models in ID (FigOU) and in 2D 
(Fig l25j) . The parameter a is taken to be a = 1. 



As can be seen from the figures, the form of the 
magnitude distribution changes only httle after 
the mainshock except that the weight of large 
events is decreased appreciably, particularly in 2D. 
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Figure 25: The time-resolved local magnitude distri- 
bution of the 2D BK model after the mainshock of 
> 1-j.c = 5. 



4 Summary and discussion 

In summary, we studied the spatiotemporal cor- 
relations of earthquakes both by the analysis of 
real earthquake catalog of Japan and by numerical 
computer simulations of the spring-block model in 
ID and 2D. Particular attention was paid to the 
magnitude distribution, the recurrence-time dis- 
tribution, the time evolution of the spatial dis- 
tribution of seismicity before and after the main- 
shock, and the time evolution of the magnitude 
distribution before and after the mainshock. Cer- 
tain eminent features of the spatiotemporal corre- 
lations, including foreshocks, aftershocks, swarms 
and doughnut-like seismic pattern, were discussed 
in some detail. 

In our numerical simulations of the BK model, 
particular attention was paid to the issue how 
the statistical properties of earthquakes depend on 
the frictional properties of earthquake faults. We 
have found that when the extent of the velocity- 
weakening property is suppressed, the system 
tends to be more critical, while, as the velocity- 
weakening property is enhanced, the system tends 
to be more off-critical with enhanced features of 
characteristic earthquakes. 

Overall, the BK model tends to exhibit more 
characteristic or off-critical statistical properties, 
particularly for large earthquakes, than the real 
seismicity which exhibits much more critical sta- 
tistical properties. This discrepancy between the 



model and the real seismicity has been recognized 
for some time now, but its true cause has remained 
to be unclear. 

First, we need to recognize that the earthquake 
catalog is taken not for a single fault, but over 
many faults. There exists a suspicion that, even 
if the property of a single individual fault is more 
or less characteristic the property obtained after 
averaging over many faults each of which has dif- 
ferent characteristics, becomes apparently char- 
acteristic as a whole. If this is really the case, 
the real observation is not necessarily inconsistent 
with the observation for the BK model, since the 
BK model deals with the property of a single uni- 
form fault. There is a claim that the extent of 
the criticality of earthquakes might depend on the 
type of earthquake faults, i.e., a matured fault 
with relatively regular fault zone behaves more 
characteristic, while an inmatured fault with rel- 
atively irregular fault zone behaves more critical 
52 . Difficulty in testing such a hypothesis is that 
the statistical accuracy of events for a single fault 
is rather limited, particularly for large events. We 
need to be careful because, when the event num- 
ber is not sufficient, an apparent deviation from 
the criticality might well arise simply due to the 
insufficient statistics, pretending a characteristic 
earthquake. 

Other possibility is that, since smaller earth- 
quakes are more or less critical even in the BK 
model, in real seismicity, the critical behavior 
might be limited to moderately large earthquakes 
which are contained in enough number in the 
earthquake catalog, while very large earthquakes, 
which are very few in number in the catalog, 
might be more or less characteristic. Anyway, the 
question of either critical or characteristic earth- 
quakes is one of major fundamental questions left 
in earthquake studies. 

The BK model was found to exhibit several 
intriguing precursory phenomena associated with 
large events: Preceding the mainshock, the fre- 
quency of smaller events is gradually enhanced, 
whereas just before the mainshock it is suppressed 
only in a close vicinity of the epicenter of the up- 
coming mainshock (the Mogi doughnut). The ap- 
parent B-value of the magnitude distribution in- 
creases significantly preceding the mainshock. On 
the other hand, the Omori law of aftershocks is 
not observed in the BK model. 

Some of these precursory phenomena observed 
in the BK model are also observed in real earth- 
quake catalog, but some of them are not. For ex- 
ample, the enhancement of foreshock activity is 
observed in common both in the BK model and in 
the JUNEC catalog. By contrast, the doughnut- 
like quiescence generically observed in the BK 



model is not observed in standard earthquakes in 
the JUNEC catalog, although it is observed in cer- 
tain earthquake swarms. 

Here, in order to make a further link between 
the BK model and the real world, it might be of 
some interest to estimate various time and length 
scales involved in the BK model. For this, we 
need to estimate the units of time and length of 
the BK model in terms of real- world earthquakes. 
Concerning the time unit uj~^, we estimate it via 
the rise time of large earthquakes, ~ tt/o;, which is 
typically about 10 seconds. This gives an estimate 
of uj~^ ~ 3 sec. Concerning the length unit L, we 
estimate it making use of the fact that the typi- 
cal displacement in large events of our simulation 
is of order one L unit, which in real- world large 
earthquakes is typically 5 meters. Then, we get 
L ~ 5 meters. Since the loading rate i^' associated 
with the real plate motion is typically 5 cm/year, 
the dimensionless loading rate = N/{Llo) is es- 
timated to be ~ 10^^. If we remember that the 
typical mean recurrence time of large events in 
our simulation is about one unit of the mean 
recurrence time of our simulation corresponds to 
100 years in real world. 

In our simulation of the BK model, the 
doughnut-like quiescence was observed before the 
mainshock at the time scale of, say, tv "S 10~^. 
This time scale corresponds to about 1 year. In 
our simulation, the doughnut-like quiescence was 
observed in the region only within a few blocks 
from the epicenter of the mainshock. To give the 
corresponding real- world estimate, we need the 
real-world estimate of our block size a. In the BK 
model, the length scale a is entirely independent 
of the length scale L, and has to be determined 
independently. We estimate a via the typical ve- 
locity of the rupture propagation, law, which is 
about 3 km/sec in real earthquakes. From this 
relation, we get a ~ 3 km. The length scale as- 
sociated with the doughnut-like quiescence is then 
estimated to be 3~6 km in radius. If we remember 
that the rupture size of large events in our simula- 
tion with = fic = 3 was about 60 blocks, the size 
of the rupture size of large earthquakes of our sim- 
ulation corresponds to 180 km and more. This is 
comparable to the size of the rupture zone of real 
earthquakes of their magnitude eight. Hence, the 
large events in the BK model might correspond to 
exceptionally large earthquakes in real seismicity, 
which might explain the reason, at least partially, 
why the deviation from the GR law observed in the 
BK model at larger magnitudes is hardly observed 
in real seismicity. In real faults, the possible max- 
imum size of the rupture zone might be limited 
by the fault geometry, i.e., by the boundary of a 



fault. 

Of course, it is not a trivial matter at all how 
faithfully the statistical properties as observed for 
the BK model represent those of real earthquakes. 
We should be careful not to put too much meaning 
to the quantitative estimates given above. 

As an other precursory effect, the change of the 
apparent 5-value is observed in the BK model, 
i.e., an increase in ID or an initial decrease fol- 
lowed by a subsequent increase in 2D. In the 
JUNEC catalog, the i?-value turns out to decrease 
prior to the mainshock consistently with several 
other observations for real faults. Meanwhile, an 
increase of the i?-value, which is similar to the 
one observed in the ID BK model, was reported 
in some real faults. Thus, to elucidate the detailed 
mechanism behind the change of the i?- value pre- 
ceding the large event is an interesting open ques- 
tion. 

One thing seems to be clear: Much needs to 
be done before we understand the true nature of 
earthquakes. We hope that further progress in 
statistical-physical approach to earthquakes, com- 
bined with the ones in other types of approaches 
from various branches of science, would eventually 
promote our fuller understanding of earthquakes. 

The author is thankful to Mr. T. Mori and Mr. 
A. Ohmura for their collaboration and discussion. 
He is also thankful to Prof. B. Chakrabarti for or- 
ganizing the MOE conference and for giving me an 
opportunity of presenting a talk there and writing 
this article. 



References 



C.H. Scholz, Nature 391, 3411 (1998). 

C.H. Scholz, The Mechanics of Earthquakes 
and Faulting (Cambridge Univ. Press, 1990). 

P. Bak, C. Tang and K. Wiesenfeld, Phys. 
Rev, Lett. 59, 381 (1987). 

P. Bak and C. Tang, J. Geophys. Res. 94, 
15635 (1989). 

R. Burridge and L. Knopoff, Bull. Seismol. 
Soc. Am. 57 (1967) 3411. 

T. Mori and H. Kawamura, Phys. Rev. Lett. 
94, 058501 (2005). 

T. Mori and H. Kawamura , J. Geophys. Res., 
m press (|physics/050"4218l). 



T. Mori and H. Kawamura, unpublished. 
H. Kawamura and T. Mori, unpublished. 



[10] J.M. Carlson, J.S. Langer, B.E. Shaw and C. 
Tang, Phys. Rev. A44, 884 (1991). 

[11] J.M. Carlson, J.S. Langer and B.E. Shaw, 
Rev. Mod. Phys. 66, 657 (1994). 

[12] B.E. Shaw, J. Geophys. Res. 100, 18239 
(1995). 

[13] C.R. Myers, B.E. Shaw and J.S. Langer, 
Phys. Rev, Lett. 77, 972 (1996). 

[14] J. Dietrich, J. Geophys. Res. 84, 2161 (1979). 

[15] A. Ruina, J. Geophys. Res. 88, 10359 (1983). 

[16] S.T. Tse and J.R. Rice, J. Geophys. Res. 91, 
9452 (1986). 

[17] N. Kato, J. Geophys. Res. 109, B12306 
(2004). 

[18] A. Ohmura and H. Kawamura, unpublished. 

[19] J.M. Carlson, Phys. Rev. A44, 6226 (1991). 

[20] J.M. Carlson and J.S. Langer, Phys. Rev. 
Lett. 62, 2632 (1989); Phys. Rev. A40, 6470 
(1989). 

[21] J.M. Carlson, J. Geophys. Res. 96, 4255 
(1991). 

[22] J. Schmittbuhl, J.-P. Vilotte and S. Roux, J. 
Geophys. Res. 101, 27741 (1996). 

[23] B.E. Shaw, J.M. Carlson and J.S. Langer, J. 
Geophys. Res. 97, 479 (1992). 

[24] C.R. Myers and J.S. Langer Phys. Rev. E47, 
3048 (1993). 

[25] B.E. Shaw, Geophys. Res. Lett. 21, 1983 
(1994). 

[26] R. De and G. Ananthakrisna, Europhys. Lett. 
66, 715 (2004). 

[27] J. Xia, H. Gould, W. Klein and J.B. Run- 
die, Phys. Rev. Lett. 95, 248501 (2005); 
cond-mat/0601679, 

[28] H. Nakanishi, Phys. Rev, A41, 7086 (1990). 

[29] K. Ito and M. Matsuzaki, J. Geophys. Res. 
95, 6853 (1990). 

[30] S.R. Brown, C.H. Scholz and J.B. Rundle, 
Geophys. Res. Lett. 18, 215 (1991). 

[31] Z. Olami, H.J. Feder and K. Christensen, 
Phys. Rev, Lett. 68, 1244 (1992). 



[32] S. Hergarten and H. Neugebauer, Phys. Rev. 
E61, 2382 (2000). 

[33] S. Hainzl, G. Zoller and J. Kurths, J. Geo- 
phys. Res. 104, 7243 (1999); Geophys. Res. 
Lett. 27, 597 (2000). 

[34] A. Helmstetter, S. Hergarten and D. Sor- 
nette, Phys. Rev, Lett. 88, 238501 (2002); 
Phys. Rev. E70, 046120 (2004). 

[35] H. Tsuruoka, Abstract of the annual meet- 
ing of the Seismological Society of Japan, p04 
(1997). 

[36] G.L. Vasconcelos, Phys. Rev. Lett. 76, 4865 
(1996). 

[37] M. S. Vieira, G.L. Vasconcelos and S.R. 
Nagel, Phys. Rev. E47, R2221 (1993). 

[38] I. Clancy and D. Corcoran, Phys. Rev. E71, 
046124 (2005). 

[39] T.C. Hanks, J. Geophys. Res. 84, 2235 

(1979) . 

[40] D.J. Andrews, J. Geophys. Res. 85, 3867 

(1980) . 

[41] K. Aki, in Earthquake Prediction, An In- 
ternational Review, (ed. D.W. Simpson and 
P.G. Richards, American Geophysical Union, 
Washington D.C.) p.566 (1981). 

[42] B.K. Chakarabarti and R.B. Stinchcombe, 
Physica, A270, 27 (1990); P. Bhat- 
tacharyya, Physica, A348, 199 (2005); P. 

Bhattacharyya, A. Chatterjee and B.K. 

Chakarabarti, [physics/0 510038 , 

[43] S.P. Nishenko and R. Buland, Bull. Seismol. 
Soc. Am. 77, 1382 (1987). 

[44] P. Bak, K. Christensen, L. Danon and T. 
Scanlon, Phys. Rev. Lett. 88, 178501 (2002). 

[45] A. Corral, Phys. Rev. Lett. 92, 108501 
(2004); Phys. Rev. E68, 035102 (2003). 

[46] K. Mogi, Bull. Earthquake Res. Inst. Univ. 
Tokyo 47, 395 (1969); Pure Appl. Geophys. 
117, 1172 (1979). 

[47] S.L. Pepke, J.M. Carlson and B.E. Shaw, J. 
Geophys. Res. 99, 6769 (1994). 

[48] V.G. Kossobokov and J.M. Carlson, J. Geo- 
phys. Res. 100, 6431 (1995). 

[49] S. Suehiro, T. Asada and M. Ohtake, Papers 
Meteorol. Geophys. 15, 71 (1964). 



[50] S.C. Jaume and L.R. Sykes, Pure Appl. Geo- 
phys. 155, 279 (1999). 

[51] W.D. Smith, Nature 289, 136 (1981). 

[52] S.G. Wesnousky, Nature 335, 340 (1988); 
M.W. Stirling, S.G. Wesnousky and K. Shi- 
mazaki, J. Goephys. Res. 124, 833 (1996). 



